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Abstract: A novel approach is presented for computing optode placements 
that are adapted to specific geometries and tissue characteristics, e.g., in op- 
tical tomography and photodynamic cancer therapy. The method is based on 
optimal control techniques together with a sparsity-promoting penalty that 
favors pointwise solutions, yielding both locations and magnitudes of light 
sources. In contrast to current discrete approaches, the need for specifying 
an initial set of candidate configurations as well as the exponential increase 
in complexity with the number of optodes are avoided. This is demonstrated 
with computational examples from photodynamic therapy. 
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1. Introduction 

In biological and medical sciences, diagnostic and therapeutic instrumentation based on visible, 
near-infrared or near-ultraviolet light (referred to as optical imaging) is of special interest as it 
is often non-invasive, the cost for the equipment is moderate, and data acquisition is usually fast 
(compared to, e.g., MRI). However, the high scattering coefficient of biological tissues in the 
visible spectrum can be a limiting factor for this technique as it causes photons to propagate in 
a non-deterministic manner. This in turn makes measurements of superficial structures and the 
selective illumination of deeper regions significantly more complicated. Due to the stochastic 
nature of the photon paths, the optimal placement of the optodes is not trivial except for simple 
regular geometries such as cylinders or spheres. 

In recent years, there has thus been increased interest in computation-driven optimization of 
optical hardware operating in strongly scattering tissues. For example, [1] used a singular value 
analysis (SVA) of the sensitivity matrix (there called "weight matrix") to compare measurement 
setups which differed either in the optode spacing or the measurement type (reflectance vs. 
transmittance setup). In [2], two optode configurations for a hybrid MRI/DOT measurement 
device were compared, and the number of singular values above a certain threshold was used as 
a quality criterion. The work in [3] applied the SVA approach to assess different fluorescence 
tomography setups in two dimensions, while [4] performed comparisons of three-dimensional 
setups. 

A similar problem is faced in photodynamic therapy (PDT) [5], which is used for dermatolog- 
ical and oncological treatments (e.g., esophageal cancers, especially at a stage where surgical 
intervention is not indicated). It appears attractive to extend PDT to other carcinomas on epi- or 
endothelial surfaces, and clinical trials are earned out for, e.g., cervix carcinomas. Some other 
potential candidates, like mesotheliomas of the thoracic cavity (which are typically difficult 
to treat), represent a special challenge. A major problem is the design of an appropriate light 
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applicator. In the esophageal cavity, which has a simple geometry, the laser light can be read- 
ily applied using a cylindrical scattering device. In contrast, the geometry of the intrathoracic 
cavity is very complex. The application of a standard esophageal applicator is not recommend- 
able because its area of treatment is small and curved and therefore the illumination would be 
inhomogeneous. Homogeneity of the irradiation is, however, a crucial design criterion, as in- 
homogeneities can lead to locally ineffective treatment on the one hand and local overdoses on 
the other hand. This can have serious consequences, including lethal overdoses [6]. 

In the past, some designs for flexible light diffusers have been proposed, where a typical 
solution is to use cylindrical or spherical diffusers in a bag filled with a scattering medium. 
These are applicable after pneumonectomia, if blood accumulations at the surface of the bags 
are prevented by continuous rinsing [7-10]. Some regions like the sinus diaphragmaticus are 
difficult to access and have to be illuminated separately. This can be achieved with wedge- 
shaped illuminators [11]; however, positioning of these illuminators and homogenization of the 
illumination is difficult. Another approach is to fill the thorax with a biologically non-hazardous 
scattering medium (e.g., with intralipid), which can also be used for rinsing to avoid blood 
accumulation. Typically, a spherical diffuser is used for illumination. However, it is difficult to 
control the dose rate with this approach. This method has also been applied to cases where lung 
tissue was not resected [8]. In general, both methods try to achieve homogeneous fluence using 
real-time dosimetry and manual repositioning of the light diffuser. An interesting alternative 
consists in textile-based diffusers, where special optical fibers are integrated in a textile. These 
diffusers are very flexible but suffer from inhomogeneous illumination and a low transmission 
rate [12, 13]. Recently so called "light blankets" with arrays of cylindrical diffusers [14] or a 
spirally-wound side-glowing fiber [15] embedded in a bag filled with intralipid were presented. 
They are easy to fabricate but still show inhomogeneities, especially at the corners. Due to the 
need for a homogeneous fluence rate, it is of great interest to optimize the placement of these 
fibers to obtain improved illumination using minimal energy. 

The purpose of this work is thus to present a general approach to compute adapted optode 
locations for different geometries, tissue types, and applications. The method is based on con- 
sidering this task as an optimal control problem for a partial differential equation describing 
the diffusion of photons in a strongly scattering medium, where the location of optodes are 
modeled as a continuous "source field". The crucial step — first proposed in [16] — is to include 
a penalty term that favors pointwise solutions. In this way, both locations and magnitudes of 
the light sources to be placed are obtained in a single step. The main advantage of this ap- 
proach over previously published — discrete — methods is that no initial maximal or minimal 
configuration needs to be specified (although an allowable region can be enforced), and that 
a combinatorial problem with exponential complexity is avoided. In addition, the algorithm is 
not based on stochastic (e.g., Monte Carlo) methods but is fully deterministic, which eases the 
verification of the outcome significantly. Finally, the approach is flexible and can incorporate a 
wide variety of objective criteria (e.g., photon flux over a given boundary section) by changing 
the target functional. The proposed approach is demonstrated in the context of optimizing the 
illumination pattern in the photodynamic treatment of intrathoracic cancer. 

The rest of this work is organized as follows. In section 2.1, the specific mathematical model 
for photon diffusion is given. Section 2.2 presents our optimal control framework for optode 
placement, whose numerical solution is described in section 2.3. The setup of the numerical 
experiments used for demonstrating our approach can be found in section 3, and the results are 
presented in section 4. A discussion of the proposed method in section 5 concludes the work. 
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2. Theory 

During photodynamic treatment of cancer, a photosensitizer such as Photofrin is injected in- 
travenously. Afterwards, the cancerogeneous site is illuminated with red to near-infrared light 
from sources in a diffuser which is applied directly on the region of interest, i.e., in the in- 
trathoracic cavity. The absorption of energy by the photo-activable drug leads to the formation 
of cytotoxic singlet oxygen, which destroys cancer cells selectively. The challenge is to homog- 
enize the light intensity as both under- and overexposure can lead to ineffective treatment [17]. 

2.1. Mathematical model 

We use the diffusion approximation of the radiative transfer equation to model the steady state 
of light propagation in a scattering medium [18]. This leads to a stationary elliptic partial dif- 
ferential equation for the photon distribution ^ G //' (H), 

-V ■{K{x)V(p{x))+Ha{x)(p{x)=q{x) inn, ^^^^ 
K{x)n(x) ■V^[x) + p^{x) =Q onT. 

The geometry of the object is given by the domain Q. C W', d £ {2,3} being the number 
of spatial dimensions, with boundary Y whose outward normal vector is denoted by n. The 
medium is characterized by the absorption coefficient pLa, the reduced scattering coefficient /i,, 

and the diffusion coefficient K ~ (M„ + jU,')] ■ The coefficient p models the reflection of a 
part of the photons at the boundary due to a mismatch in the index of refraction. Finally, the 
source term q models the light emission of the embedded optodes. 

For the optimal control approach, we also require the solution p e //' (i2) of the adjoint 
equation 

-y-{K{x)yp{x))+pLa{x)p{x)^f{x) inn, ^^^^ 
K{x)n{x)-Vp{x)+pp{x)=Q onF 

for given / G L^(n). Both equations should be understood in the weak sense. 

2.2. Optode placement optimization 

Since optodes act as discrete light sources, the source term can be modeled as q{x) ~ 
Ylj=\ qjS{x — Xj) for qj G IR+ and xj G H, 1 < y < A^, where 8 denotes the Dirac distribution 
(i.e., / fd8{x) — /(O) for all continuous functions /). A straightforward approach for optimiz- 
ing the placement of the optodes (as was done, e.g., in [19]) would identify a set of M ^ 
possible optode locations xi ,.. .,xm and chose the best locations such that a certain perfor- 
mance criterion J{q) is minimized. The corresponding optimal source magnitudes qj would 
then be computed in a second step. 

To avoid the combinatorial complexity of this discrete approach, instead of specifying the 
optode locations beforehand, we optimize the (distributed) source term q directly while adding 
a penalty term that promotes sparsity of q, i.e., smallness of its support {x G : q{x) ^ 0}. 
This has the added advantage that the number of optodes need not be specified in advance. 
Since we are looking for point sources, this requires searching for q in the space of regular 
Borel measures (which includes the Dirac distribution). Following [20], we are thus led to the 
optimization problem 

min y(^) + a||^||.^, 
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where (il) is the space of regular Borel measures, i.e., the dual of the space Co{Q.) of con- 
tinuous functions with compact support on £2, with norm 

1I?IU= sup / fdq, 
ll/llc<l 

which reduces to 

= / \q{x)\dx= ll^ll^i 
Jn 

for q G L'(£2). This is related to the well-known fact that L' norms promote sparsity in op- 
timization. The penalty parameter a controls the sparsity of the solution: The larger a, the 
smaller the support of q. 

Motivated by the application in PDT, we chose as performance criterion the deviation from 
a constant illumination z in an observation region (Oo CQ. such that y(^) := 4||<p|(»„ — z||?2/,, \> 
where (p\(Og denotes the restriction of (p to cOo- Due to the linearity of the forward problem, 
we can take z = 1 Wm^^ without loss of generality. After optimization, the magnitude of the 
resultant sources can be linearly scaled to achieve the required illumination z- In addition, we 
restrict the possible light source locations to a control region (Og C £1, which does not overlap 
with the observation region a)„ (i.e., a)^, n a)„ = 0), and enforce non-negativity of the source 
term q (which represents the optodes). This leads to the following optimization problem: 

1 T 

, min 7^\\(P\o>„-z\\l2(i^) + (x\\q\\.jr{m„) subject to (2.1) and ^ > 0. (2.3) 

<peH^{Q.),qe^{o)g) ^ ^ 

It was shown in [21] that this problem has a solution cf G ^(o)^), which can be approximated 
by a sequence of functions qy^lJia^) for 7^ oo satisfying 

qy + 7min(0, + a) = 0, (2.4) 

where py is the solution of (2.2) with right hand side f :^ (py — z and (py is the solution of (2.1) 
with right hand side q^. Equation (2.4) can be solved using a semismooth Newton method which 
is superlinearly convergent; see [21]. To globalize the Newton method and closely approximate 
the solution q* of (2.3), we use a continuation scheme in 7 where we iteratively solve the 
problem for an increasing sequence 7„, using the previous solution as initial guess. 

2.3. Finite element discretization 

The discretization needs to account for the fact that the functions qy converge to measures as 
7 increases. We therefore employ the finite element discretization proposed in [22], where the 
photon density <py and the adjoint variable py are discretized using piecewise linear elements 
on a given triangulation T, while the source term qy is discretized using linear combinations of 
Dirac distributions centered at the interior nodes jc,, 1 < / < N{T), of T: 

N{T) 

!=1 

In practice, the number of nodes N{T) will be determined by the need to resolve the geometry 
of the domain and the required accuracy of the solution of the forward model (2.1). Although 
further refinement of the triangulation increases the number of possible optode locations, the 
sparsity-promoting property of the minimized functional discourages placing additional op- 
todes. In fact, it was shown in [22] that for a given discretization of the forward model, the 
computed sources (for 7 ^ 00) are optimal among all (non-discretized) measures. 
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Algorithm 1 Semismooth Newton method with continuation 



1: for wj = 1, . . . ,m* do 

2: set7=2('"-'), = = = 0 

3: foryt = 0,...,/t* do 

4: solve (2.5) for 77'^+' 

5: compute from (2.6) 

6: if£/*+i =£/*^then 

7: set^('") = 7min(0,p'-+'|fi)^ + a) 

8: break 

9: end if 

10: end for 

11: end for 



Since the linear finite element basis functions form a nodal basis, the right hand side in the 
weak formulation of (2.1) for a piecewise linear basis function Cj becomes 

N(T] 

i^Y^^j) = L qi{S{x-Xi),ej) ^qj, 

!=1 

i.e., the mass matrix is the identity. Introducing the stiffness matrix A corresponding to (2. 1) and 
the observation mass matrix M,, with entries M,y = ejCjcix, we obtain the discrete optimality 
system 

A(py — qy = Q, 

^ qy + Ymin{0,py\a^^ + a) = 0, 

Eliminating qy using the last equation and applying a semismooth Newton method, cf. [21], we 
have to solve for (^*^+' jp'^^' ) the block system 



A DA f(p''+^\ _ f-ad'' 



-Ma A J V^^"~ 

where is a diagonal matrix with the entries of the vector d'' 



(2.5) 



if(/k).<-«> (2.6) 
^ [0 else, 

on the diagonal. It can be shown that the semismooth Newton method has converged once 
of'^+i = d'^ holds. After the final has been computed, the corresponding control can be ob- 
tained from (2.4). The complete procedure is given in Algorithm 1. 

3. Materials and methods 

The optimization algorithm described in section 2.2 is implemented in Python using the open 
source finite element library FEniCS [23]. The parameters in Algorithm 1 are set to m* — 34 
(such that 7* « 10'") and k* — 20. To model a textile-based diffuser, the material parameters 
in (2.1) are taken as jXa — 10^^ mm^^ /i' = 10^^ mm^^ and p = 0.1992. The influence of the 
parameter a is illustrated by comparing the results for different values of a specified below. 

The meshes for the light diffusers containing the optodes are created with the commercial 
mesh generator Hypermesh'"^. To demonstrate the behavior of the optimization algorithm for 
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(b) Double-curved models. 



Fig. 1. Two-dimensional model geometries (numbers denote curvature K). 



different geometries, we first consider simple two-dimensional spline models which represent 
the cross-section of an infinitely long pad. This geometry mimics that of an array of parallel 
cylindrical diffusers embedded in a scattering substrate. Five single-curved models and four 
double-curved models with increasing curvature K were created as shown in Fig. 1. The di- 
mensions correspond approximately to a width of 10 mm and height of 120 mm. In all cases, 
the region a)„ in which the illumination should be homogenized are the left and right outer 
lines (indicated in orange in Fig. 1). The region o^y where optodes are allowed to be placed is 
a single line equidistant from both (indicated by a dashed line in Fig. 1). The meshes for the 
single-curved models of curvature K = 5, 10, 20, 40, and 60 consist of 61038, 61789, 67160, 
80664, and 105322 finite elements, respectively. The double-curved models of curvature K" = 5, 
10, 15, and 20 are comprised of 62349, 70735, 821 19, and 104220 finite elements, respectively. 

The photodynamic treatment is simulated by embedding the light diffuser model in the in- 
trapleural space of a realistic three-dimensional human thorax model that is constructed from 
a stack of CT images. The approximate dimensions are: height 100 mm, width 150mm, thick- 
ness 10 mm. The observation region (Oo is defined as the outer and inner surface of the model, 
and cOii is an interior manifold equidistant from both (see Fig. 2; cOg is indicated in purple). The 
generated mesh consists of 81770 elements. 

The results are evaluated quantitatively for different values of the sparsity-controlling pa- 
rameter a. The coefficient of variation c,, of the resultant photon density (py over the ob- 
servation region a)„ and the number of sources after the optimization procedure serve as 
quality measures. For the latter, the nodes in the control region (Og satisfying > 10^'^ are 
counted. We compare the results for a £ {0.1, 0.01 , 0.001 } for the two-dimensional models and 
a e {0.2,0.4, . . . , 1.8} for the three-dimensional model. 

4. Results 

The quantitative results for the two-dimensional geometries are given in Table 1 for the single- 
curved models and in Table 2 for the double-curved models. As can be seen by comparing the 
number of active nodes with the total number of nodes for each model, the algorithm in- 
deed produces discrete sources that can be used as optode positions. The obtained coefficients 
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(a) Top oblique view 




(b) Front oblique view (c) Side oblique 

view 



Fig. 2. Three-dimensional model. The admissible manifold ft)^ for optodes is indicated in 
purple. 

of variation c,, indicate that a homogeneous illumination of the desired region is possible at 
least for a < 0.1, demonstrating the feasibility of the proposed approach. The robustness of 
the algorithm with respect to geometry is illustrated by the fact that the achieved variations do 
not depend very much on the curvature. It can also be observed how the penalty parameter a 
determines the tradeoff between the number of active optodes and the homogeneity of the illu- 
mination in the region of interest: larger values of a yield fewer optodes but less homogeneous 
illumination, again independent of curvature. 

The qualitative behavior of the computed sources for each value of a is shown in Fig. 3(a) 
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and Fig. 3(b) for a representative single-curved (fc = 20) and double-curved model (k = 15), re- 
spectively, where the relative strength of the sources is coded by height. (Note when comparing 
Tables 1 and 2 with Fig. 3 that neighboring active nodes appear as a single peak and thus can 
be taken as a single optode.) While for the single-curved model and a = 0.1, the distribution of 
optodes agrees well with the intuitive choice of equally spaced optodes of approximately equal 
magnitude, the other values indicate that a better illumination can be achieved with stronger 
sources towards the tips of the model. It should be pointed out that even in the former case, the 
number of optodes to be distributed is not obvious. For the double-curved models, the results 
indicate that optodes should be placed preferentially in regions where the curvature changes. 

For reference. Fig. 4 shows the corresponding photon densities (py (in Wm^^, normalized 
to unit mean) plotted along part of the observation region (left line in Fig. 1), illustrating how 
the parameter a and the model geometry influence the homogeneity of the illumination in this 
region. As expected, photon fluence shows the most pronounced inhomogeneities close to the 
borders. In the case of the single curved model, a nearly sinusoidal ripple pattern arises in more 
than 80 % of the target region, while in the double curved model the ripple is superimposed on 
a step-profile with the steps located approximately at the zero-crossing points of the curvature. 
With a = 0.1, the peak-peak fluctuations are still around 40 % of the mean value even far away 
from the borders, which may be considered as unsatisfactory. However, when decreasing alpha 
to 0.01 or less, the ripple remains within a few percent, which is sufficient, especially when 
comparing this value to other sources of fluctuations of the irradiation such as local absorp- 
tion changes by tissue inhomogeneities, bleeding, or inhomogeneities of the distribution of the 
photosensitizer. 

The quantitative results for the three-dimensional model are shown in Table 3. For a = 
1.8, no controls are placed and thus the photon density is zero. This is consistent with the 
theory, which predicts that there is a threshold value for a above which the optimal control is 



Table 1. Results for single-curved models. Shown are the number N of active nodes and 
the coefficient of variation c,, of the photon density in the observation domain for different 
curvatures K and values of a. 
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Table 2. Results for double-curved models. Shown are the number A' of active nodes and 
the coefficient of variation c,, of the photon density in the observation domain for different 
curvatures K and values of a. 
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1.51e-2 


2.03e-l 


2.88e-2 


2.45e-2 


2.24e-l 


3.27e-2 


2.94e-2 


4.92e-l 


3.47e-2 


3.09e-2 



Table 3. Results for three-dimensional model. Shown are the number A' of active nodes and 
the coefficient of variation c,. of the photon density in the observation domain for different 
values of a. 



a 1.8 


1.6 


1.4 


1.2 


1.0 


0.8 


0.6 


0.4 


0.2 


N 0 


12 


150 


250 


333 


409 


498 


637 


884 


Cv — 


1.856-1-0 


5.64e-l 


3.59e-l 


2.65e-l 


2.04e-l 


1.56e-l 


1.13e-l 


6.72e-2 
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(a) Single-curved model (K = 20) (b) Double-curved model {K= 15) 



Fig. 4. Photon densities (py (in Wm^^, normalized to unit mean) plotted along part of 
the observation region (left line in Fig. 1) for representative single-curved and double- 
curved models for three different values of alpha (from top to bottom: a = 0.1, a = 0.01, 
a = 0.001). 



identically zero; cf. [22, Proposition 2.2]. 

Figure 5 shows location and magnitude (color coded) of the computed optodes and the cor- 
responding photon densities (in Wm^^, normalized to unit mean) for a = 1.2, a = 0.8, and 
a ~ 0.4. Due to the nonuniform curvature of the model, a homogeneous illumination is harder 
to achieve than in the two-dimensional case, especially at the borders of the target region. How- 
ever, for a < 1.2, the inhomogeneities in the interior are usually within 10 %, and the few hot 
spots of 30 % would still be acceptable. Although of course the specific placement may be dif- 
ficult to realize in practice, the qualitative distribution can be useful information in the initial 
design process. 




(c) optodes, a = 0.8 (d) photon density, a = 0.8 




(e) optodes, a = 0.4 (f) photon density, a = 0.4 

Fig. 5. Optode positions and magnitudes (left) and photon densities (right; in Wm^^, nor- 
malized to unit mean) for the three-dimensional model and three different values of a. 

5. Discussion 

The proposed approach is able to generate reasonable optode configurations adapted to specific 
geometries, even in situations (such as complex three-dimensional models) where optimal se- 
tups are not intuitively obvious. Our method also yields relative strengths of the optodes to be 
placed, which would otherwise have to be computed in a separate step. Fuithermore, the algo- 
rithm is deterministic and does not require a-priori knowledge such as an initial set of candidate 
locations or the number of optodes after optimization, which on the contrary is provided by our 
approach. The method can be used as a tool during the initial design process to estimate the 
number of sources required as well as their location and relative strengths. 

By formulating the optode placement problem as a continuous optimization problem, the 
combinatorial complexity inherent in discrete approaches is avoided. This is critical for achiev- 



ing an efficient optimization technique and — to our knowledge — has not been presented before 
in the context of diffuse optical imaging. As an example, our Python implementation required 
about three minutes on a MacBook Pro (2. 16 GHz Intel Core2 Duo with 2 GByte RAM) for the 
single-curved model with K" = 5. Our approach could therefore also be used in an interactive 
setting, where the engineer will adapt design parameters, such as the optical coefficients of the 
diffuser, based on the outcome of an optimization run. 

While the number of desired optodes is correlated with the penalty parameter a, it is not 
directly controllable. This drawback is analogous to the problem of finding the "best" regular- 
ization parameter in image reconstruction (e.g., for diffuse optical tomography), where typically 
the determination of the parameter is left to the user or is based on heuristics. Certainly, one 
could think about finding a good parameter through successive optimization runs, e.g., with 
decreasing values of a if the user specified an upper bound on the number of optodes. 

The achieved results are satisfactory from the mathematical point of view; but of course they 
should also be discussed in an engineering context. In particular, it may be difficult to place 
many sources in a more or less irregular pattern. The number of optodes depends on the required 
uniformity of the surface fluence. A reasonable value in practice would be a CV of 0.05. Table 
1 shows CVs (for single-curved pads) below 0.03 for around 50 optodes, but up to 0.25 for less 
than 25 optodes. The numbers for the double-curved pad are only slightly greater This means 
that the between 40 and 50 required sources can be expected. In practice, such a design can 
be approximated comparatively easily with parallely arranged cylindrical polymer diffusers of 
sufficiently small radius, which are fed by individual optical fibers. Instead of a fixed grid of 
diffusers, one can imagine dense bundles of uniformly spaced diffusers where only those close 
to the optimal positions are connected to the laser source. This would allow a very flexible use 
and adaptive homogenization of the fluence dependent on the individual anatomical situation 
(e.g., curvature), which is certainly desirable in the context of a personalized optimization. Such 
a concept can be realized by using fiberoptic switches with many channels. In three dimensions, 
the sources may be fiber-coupled spherical diffusers or simply open-ended fibers. Due to the 
higher number of potential positions (637 for a CV of 0.11, see Tab. 3), the construction of 
a flexible structure here may be difficult, and pre-fabricated pads that are adapted to a certain 
anatomical target geometry appear more realistic than a truly adaptive system. 

Although a rigorous sensitivity analysis has not yet been carried out, our experience indicates 
that the computed photon density distributions are relatively robust to small perturbations of 
the optode locations and magnitudes. Similarly, we did not observe significant changes in the 
results due to small random perturbations of the optical parameters. This can be attributed to 
the linearity and the strong diffusivity of the model (2. 1). Such robustness is very important for 
practical implementations because it means that the result is not very sensitive to manufacturing 
tolerances. 

One of the main advantages of the optimal control approach is its flexibility. For example, 
it is straightforward to extend the underlying model to include, e.g., inhomogeneous material 
properties or to replace the diffusion approximation by a more complicated model such as the 
radiative transfer equation. It is also possible to consider different objective criteria such as the 
photon flux through a given (boundary or internal) surface by changing the functional J{q). In 
principle, the approach can be applied to the problem of optimal experiment design for optical 
tomography, if the objective J{q) is based on a suitable sensitivity term. However, this extension 
of our method is subject to future work. 
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